home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SVDCMP.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  7KB  |  252 lines
  1. PROCEDURE svdcmp(VAR a: glmpbynp; m,n,mp,np: integer;
  2.        VAR w: glnparray; VAR v: glnpbynp);
  3. (* Programs using routine SVDCMP must define the types
  4. TYPE
  5.    glnparray = ARRAY [1..np] OF real;
  6.    glmpbynp = ARRAY [1..mp,1..np] OF real;
  7.    glnpbynp = ARRAY [1..np,1..np] OF real;
  8. in the main routine. *)
  9. LABEL 1,2,3;
  10. CONST
  11.    nmax=100;
  12. VAR
  13.    nm,l,k,j,jj,its,i: integer;
  14.    z,y,x,scale,s,h,g,f,c,anorm: real;
  15.    rv1: ARRAY [1..nmax] OF real;
  16. FUNCTION sign(a,b: real): real;
  17.    BEGIN
  18.       IF (b >= 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
  19.    END;
  20. FUNCTION max(a,b: real): real;
  21.    BEGIN
  22.       IF (a > b) THEN max := a ELSE max := b
  23.    END;
  24. BEGIN
  25.    g := 0.0;
  26.    scale := 0.0;
  27.    anorm := 0.0;
  28.    FOR i := 1 TO n DO BEGIN
  29.       l := i+1;
  30.       rv1[i] := scale*g;
  31.       g := 0.0;
  32.       s := 0.0;
  33.       scale := 0.0;
  34.       IF (i <= m) THEN BEGIN
  35.          FOR k := i TO m DO BEGIN
  36.             scale := scale+abs(a[k,i])
  37.          END;
  38.          IF (scale <> 0.0) THEN BEGIN
  39.             FOR k := i TO m DO BEGIN
  40.                a[k,i] := a[k,i]/scale;
  41.                s := s+a[k,i]*a[k,i]
  42.             END;
  43.             f := a[i,i];
  44.             g := -sign(sqrt(s),f);
  45.             h := f*g-s;
  46.             a[i,i] := f-g;
  47.             IF (i <> n) THEN BEGIN
  48.                FOR j := l TO n DO BEGIN
  49.                   s := 0.0;
  50.                   FOR k := i TO m DO BEGIN
  51.                      s := s+a[k,i]*a[k,j]
  52.                   END;
  53.                   f := s/h;
  54.                   FOR k := i TO m DO BEGIN
  55.                      a[k,j] := a[k,j]+
  56.                         f*a[k,i]
  57.                   END
  58.                END
  59.             END;
  60.             FOR k := i TO m DO BEGIN
  61.                a[k,i] := scale*a[k,i]
  62.             END
  63.          END
  64.       END;
  65.       w[i] := scale*g;
  66.       g := 0.0;
  67.       s := 0.0;
  68.       scale := 0.0;
  69.       IF ((i <= m) AND (i <> n)) THEN BEGIN
  70.          FOR k := l TO n DO BEGIN
  71.             scale := scale+abs(a[i,k])
  72.          END;
  73.          IF (scale <> 0.0) THEN BEGIN
  74.             FOR k := l TO n DO BEGIN
  75.                a[i,k] := a[i,k]/scale;
  76.                s := s+a[i,k]*a[i,k]
  77.             END;
  78.             f := a[i,l];
  79.             g := -sign(sqrt(s),f);
  80.             h := f*g-s;
  81.             a[i,l] := f-g;
  82.             FOR k := l TO n DO BEGIN
  83.                rv1[k] := a[i,k]/h
  84.             END;
  85.             IF (i <> m) THEN BEGIN
  86.                FOR j := l TO m DO BEGIN
  87.                   s := 0.0;
  88.                   FOR k := l TO n DO BEGIN
  89.                      s := s+a[j,k]*a[i,k]
  90.                   END;
  91.                   FOR k := l TO n DO BEGIN
  92.                      a[j,k] := a[j,k]
  93.                         +s*rv1[k]
  94.                   END
  95.                END
  96.             END;
  97.             FOR k := l TO n DO BEGIN
  98.                a[i,k] := scale*a[i,k]
  99.             END
  100.          END
  101.       END;
  102.       anorm := max(anorm,(abs(w[i])+abs(rv1[i])))
  103.    END;
  104.    FOR i := n DOWNTO 1 DO BEGIN
  105.       IF (i < n) THEN BEGIN
  106.          IF (g <> 0.0) THEN BEGIN
  107.             FOR j := l TO n DO BEGIN
  108.                v[j,i] := (a[i,j]/a[i,l])/g
  109.             END;
  110.             FOR j := l TO n DO BEGIN
  111.                s := 0.0;
  112.                FOR k := l TO n DO BEGIN
  113.                   s := s+a[i,k]*v[k,j]
  114.                END;
  115.                FOR k := l TO n DO BEGIN
  116.                   v[k,j] := v[k,j]+s*v[k,i]
  117.                END
  118.             END
  119.          END;
  120.          FOR j := l TO n DO BEGIN
  121.             v[i,j] := 0.0;
  122.             v[j,i] := 0.0
  123.          END
  124.       END;
  125.       v[i,i] := 1.0;
  126.       g := rv1[i];
  127.       l := i
  128.    END;
  129.    FOR i := n DOWNTO 1 DO BEGIN
  130.       l := i+1;
  131.       g := w[i];
  132.       IF (i < n) THEN BEGIN
  133.          FOR j := l TO n DO BEGIN
  134.             a[i,j] := 0.0
  135.          END
  136.       END;
  137.       IF (g <> 0.0) THEN BEGIN
  138.          g := 1.0/g;
  139.          IF (i <> n) THEN BEGIN
  140.             FOR j := l TO n DO BEGIN
  141.                s := 0.0;
  142.                FOR k := l TO m DO BEGIN
  143.                   s := s+a[k,i]*a[k,j]
  144.                END;
  145.                f := (s/a[i,i])*g;
  146.                FOR k := i TO m DO BEGIN
  147.                   a[k,j] := a[k,j]+f*a[k,i]
  148.                END
  149.             END
  150.          END;
  151.          FOR j := i TO m DO BEGIN
  152.             a[j,i] := a[j,i]*g
  153.          END
  154.       END ELSE BEGIN
  155.          FOR j := i TO m DO BEGIN
  156.             a[j,i] := 0.0
  157.          END
  158.       END;
  159.       a[i,i] := a[i,i]+1.0
  160.    END;
  161.    FOR k := n DOWNTO 1 DO BEGIN
  162.       FOR its := 1 TO 30 DO BEGIN
  163.          FOR l := k DOWNTO 1 DO BEGIN
  164.             nm := l-1;
  165.             IF ((abs(rv1[l])+anorm) = anorm) THEN GOTO 2;
  166.             IF ((abs(w[nm])+anorm) = anorm) THEN GOTO 1
  167.          END;
  168. 1:         c := 0.0;
  169.          s := 1.0;
  170.          FOR i := l TO k DO BEGIN
  171.             f := s*rv1[i];
  172.             IF ((abs(f)+anorm) <> anorm) THEN BEGIN
  173.                g := w[i];
  174.                h := sqrt(f*f+g*g);
  175.                w[i] := h;
  176.                h := 1.0/h;
  177.                c := (g*h);
  178.                s := -(f*h);
  179.                FOR j := 1 TO m DO BEGIN
  180.                   y := a[j,nm];
  181.                   z := a[j,i];
  182.                   a[j,nm] := (y*c)+(z*s);
  183.                   a[j,i] := -(y*s)+(z*c)
  184.                END
  185.             END
  186.          END;
  187. 2:         z := w[k];
  188.          IF (l = k) THEN BEGIN
  189.             IF (z < 0.0) THEN BEGIN
  190.                w[k] := -z;
  191.                FOR j := 1 TO n DO BEGIN
  192.                v[j,k] := -v[j,k]
  193.             END
  194.          END;
  195.          GOTO 3
  196.          END;
  197.          IF (its = 30) THEN BEGIN
  198.             writeln ('no convergence in 30 SVDCMP iterations'); readln
  199.          END;
  200.          x := w[l];
  201.          nm := k-1;
  202.          y := w[nm];
  203.          g := rv1[nm];
  204.          h := rv1[k];
  205.          f := ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
  206.          g := sqrt(f*f+1.0);
  207.          f := ((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
  208.          c := 1.0;
  209.          s := 1.0;
  210.          FOR j := l TO nm DO BEGIN
  211.             i := j+1;
  212.             g := rv1[i];
  213.             y := w[i];
  214.             h := s*g;
  215.             g := c*g;
  216.             z := sqrt(f*f+h*h);
  217.             rv1[j] := z;
  218.             c := f/z;
  219.             s := h/z;
  220.             f := (x*c)+(g*s);
  221.             g := -(x*s)+(g*c);
  222.             h := y*s;
  223.             y := y*c;
  224.             FOR jj := 1 TO n DO BEGIN
  225.                x := v[jj,j];
  226.                z := v[jj,i];
  227.                v[jj,j] := (x*c)+(z*s);
  228.                v[jj,i] := -(x*s)+(z*c)
  229.             END;
  230.             z := sqrt(f*f+h*h);
  231.             w[j] := z;
  232.             IF (z <> 0.0) THEN BEGIN
  233.                z := 1.0/z;
  234.                c := f*z;
  235.                s := h*z
  236.             END;
  237.             f := (c*g)+(s*y);
  238.             x := -(s*g)+(c*y);
  239.             FOR jj := 1 TO m DO BEGIN
  240.                y := a[jj,j];
  241.                z := a[jj,i];
  242.                a[jj,j] := (y*c)+(z*s);
  243.                a[jj,i] := -(y*s)+(z*c)
  244.             END
  245.          END;
  246.          rv1[l] := 0.0;
  247.          rv1[k] := f;
  248.          w[k] := x
  249.       END;
  250. 3:   END
  251. END;
  252.